//
// nnFinEst.cc
//
// usa una rete neurale per prevedere la variazione del cambio USD / EUR
// come funzione delle variazioni dei cambi USD / EUR, GBP / EUR e JPY / EUR
// nei giorni precedenti
//
// Istruzione per l'uso:
// - caricare in un directory:
//   - nnFinEst.cc
//   - jpy-eur-1999-2005.csv  
//   - gbp-eur-1999-2005.csv
//   - usd-eur-1999-2005.csv  
// - eseguire Root:
//   - shell> root <return>
// - invocare il programma da Root chiedendo la compilazione
//   - root [1] .x nnFinEst.cc+
//

#if !defined(__CINT__) || defined(__MAKECINT__)

#include <assert.h>
#include <fstream>
#include <iostream>
#include <iomanip>
#include <string>
#include <deque>
// #include <algorithm>
// #include <functional>
#include <numeric>

#include "TSystem.h"
#include "TCanvas.h"
#include "TMath.h"
#include "TFile.h"
#include "TTree.h"
#include "TEventList.h"
#include "TMultiLayerPerceptron.h"
#include "TMLPAnalyzer.h"

#endif

//
// csvToRoot
//
// - csvFile:	file con cambio in funzione del giorno
// - rootFile:	nome file Root per i dati convertiti
// - treeName:	nome n-upla
// - nSkip:     giorni di sfasamento
// - nSkipMax:	n. massimo di giorni di sfasamento
//
// Da una sequenza temporale di cambi (es. dollaro/euro) viene generato un
// file contenente una n-upla di Root.  La serie dei cambi viene
// "regolarizzata" seguendo L.Anastasakis, N.Mort, "Applying a
// feedforward neural network for the prediction of the USD/GBP
// exchange rate", cioe':
// - si calcola logaritmo del rapporto tra:
//   - il valore del cambio al giorno i
//   - la media del valore del cambio nei 5 giorni precedenti
// - alternativamente si calcola la differenza diviso
//   la media del valore del cambio nei 5 giorni precedenti
// Questa funzione consente di generare serie di n-uple di Root
// corrispondenti a valori della variazione del cambio sfasati
// di 1,2,3...N giorni: questo facilita la costruzione di una rete
// neurale che utilizza i dati degli N giorni precedenti per predire
// il valore futuro del cambio.
//

int csvToRoot(TString csvFile, TString rootFile,
	      TString treeName = "nt",
	      UInt_t nAver=5,
	      Int_t nSkip=0, Int_t nSkipMax=0 )
{
  std::string line;
  Float_t time;
  Float_t quote;
  Float_t rwi;
  Float_t quoteDelta;
  Float_t quoteLast(0);
  // spazio per tenere i valori del cambio da mediare
  deque<Float_t> quoteBuff;

  std::ifstream input;

  input.open(csvFile);
  if (!input.is_open ()) {
    std::cerr << "Cannot open " << csvFile << endl;
  }

  //--- conta quanti valori di cambio ha il file
  int inputLines(0);
  while(std::getline(input, line)) {
    size_t firstComma = line.find(',');
    if (sscanf(line.substr(firstComma+1).c_str(), "%f", &quote) == 1) {
      inputLines++;
    }
  }
  //
  // il numero di variazioni del valore del cambio e' il numero
  // dei cambi meno uno
  //
  inputLines--;

  // std::cout << "il " << inputLines << endl;
  input.clear();
  input.seekg(0, ios::beg); 
  
  TFile out(rootFile,"RECREATE");

  int lineNum(0);
  UInt_t nConvDone(0);
  bool couldNotConvert(false);
  struct tm unixTime;
  memset(&unixTime, 0, sizeof(unixTime));

  TTree tree(treeName, treeName);
  tree.Branch("t", &time, "t/F");
  tree.Branch("q", &quoteDelta, "q/F");
  //
  // Random Walk Indicator
  // per una sequenza di valori del cambio,
  // la differenza tra l'ultimo e il primo diviso
  // il valore medio del valore assoluto delle differenze
  // rwi = (q1 - q5)/ [Sum_i abs(q(i) - q(i-1))] / (N-1)
  //
  tree.Branch("rwi", &quoteDelta, "rwi/F");

  while(std::getline(input, line)) {
    if (couldNotConvert) {
      std::cerr << "Could not convert line " << lineNum << std::endl;
      std::cerr << "-- " << line << " --" << std::endl;
    }
    lineNum++;

    size_t firstComma = line.find(',');
    if (sscanf(line.substr(firstComma+1).c_str(), "%f", &quote)
        == 1) {
      couldNotConvert = false;
      nConvDone++;

      strptime(line.substr(0,firstComma).c_str(), "%d.%m.%Y", &unixTime);
      //--- date in seconds since 1970
      time = float(mktime(&unixTime));
      //--- divide by seconds per day
      time /= 24*60*60;

      if (nConvDone > nAver) {
	quoteLast = accumulate(quoteBuff.begin(), quoteBuff.end(),0.0f)/nAver;
	// quoteDelta = 100*TMath::Log(quote/quoteLast);
	quoteDelta = 100*(quote-quoteLast)/quoteLast;

	rwi = quoteBuff[quoteBuff.size()-1] - quoteBuff[0];
	Float_t averStep;
	averStep = 0;
	for(size_t i=1; i < quoteBuff.size(); i++) {
	  averStep += abs(quoteBuff[i] - quoteBuff[i-1]);
	}
	rwi /= (averStep * sqrt(static_cast<float>(quoteBuff.size()-1)));
#if 0
	cout
	  << nConvDone-nAver << " "
	  << nSkipMax << " "
	  << nSkip << " "
	  << inputLines << " "
	  << endl;
#endif
	
	//--- riempi un numero uguale di variazioni per ogni sfasatura
	if (nConvDone > (nSkipMax - nSkip) + nAver &&
	    nConvDone < (inputLines - nSkip) + nAver) {
#if 0
      std::cout
	<< rwi
	<< " "
        << quoteBuff.size()
	<< " "
        << quoteDelta
        << " -> "
	<< quote
        << " - "
	<< time
        << " = "
        << line.substr(0,firstComma)
        << " = "
        << asctime(&unixTime)
        << std::endl;
#endif

	  tree.Fill();
	}
      }
      quoteBuff.push_back(quote);
      if (quoteBuff.size() > nAver) {
	quoteBuff.pop_front();
      }

    } else {
      //--- first line has headers
      if (lineNum>1) {
        couldNotConvert = true;
      }
    }
  }

  //--- operazioni conclusive: chiusura file di input, scrittura degli output
  input.close();
  tree.Write();
  out.Write();
  out.Close();

  return 0;
}

void nnFinEst()
{
  gROOT->Reset(); // fa il reset di tutte le variabili

  if (!gROOT->GetClass("TMultiLayerPerceptron")) {
    gSystem->Load("libMLP");
    cout
      << endl
      << "E' stata caricata la libreria ROOT per le reti neurali."
      << endl
      << endl;
  }

  //
  // converte i dati sul cambio il n-uple di Root
  // per ogni cambio vengono creati 5 file cotenenti gli stessi dati
  // ma sfasati di 0,1,2,3,4 giorni, in modo da poter utilizzare
  // all'ingresso della rete neurale i cambi dei giorni precedenti
  // per calcolare il cambio futuro
  // 

  csvToRoot("usd-eur-1999-2005.csv", "usd-eur-1999-2005.root", "usd", 5, 0, 4);
  csvToRoot("usd-eur-1999-2005.csv", "usd-eur-1999-2005-1.root", "usd1", 5, 1, 4);
  csvToRoot("usd-eur-1999-2005.csv", "usd-eur-1999-2005-2.root", "usd2", 5, 2, 4);
  csvToRoot("usd-eur-1999-2005.csv", "usd-eur-1999-2005-3.root", "usd3", 5, 3, 4);
  csvToRoot("usd-eur-1999-2005.csv", "usd-eur-1999-2005-4.root", "usd4", 5, 4, 4);

  csvToRoot("gbp-eur-1999-2005.csv", "gbp-eur-1999-2005.root", "gbp", 5, 0, 4);
  csvToRoot("gbp-eur-1999-2005.csv", "gbp-eur-1999-2005-1.root", "gbp1", 5, 1, 4);
  csvToRoot("gbp-eur-1999-2005.csv", "gbp-eur-1999-2005-2.root", "gbp2", 5, 2, 4);
  csvToRoot("gbp-eur-1999-2005.csv", "gbp-eur-1999-2005-3.root", "gbp3", 5, 3, 4);
  csvToRoot("gbp-eur-1999-2005.csv", "gbp-eur-1999-2005-4.root", "gbp4", 5, 4, 4);

  csvToRoot("jpy-eur-1999-2005.csv", "jpy-eur-1999-2005.root", "jpy", 5, 0, 4);
  csvToRoot("jpy-eur-1999-2005.csv", "jpy-eur-1999-2005-1.root", "jpy1", 5, 1, 4);
  csvToRoot("jpy-eur-1999-2005.csv", "jpy-eur-1999-2005-2.root", "jpy2", 5, 2, 4);
  csvToRoot("jpy-eur-1999-2005.csv", "jpy-eur-1999-2005-3.root", "jpy3", 5, 3, 4);
  csvToRoot("jpy-eur-1999-2005.csv", "jpy-eur-1999-2005-4.root", "jpy4", 5, 4, 4);

  TFile file("usd-eur-1999-2005.root");
  TTree* tree = (TTree*)file.Get("usd");

  tree->AddFriend("usd1", "usd-eur-1999-2005-1.root");
  tree->AddFriend("usd2", "usd-eur-1999-2005-2.root");
  tree->AddFriend("usd3", "usd-eur-1999-2005-3.root");
  tree->AddFriend("usd4", "usd-eur-1999-2005-4.root");

  tree->AddFriend("gbp", "gbp-eur-1999-2005.root");
  tree->AddFriend("gbp1", "gbp-eur-1999-2005-1.root");
  tree->AddFriend("gbp2", "gbp-eur-1999-2005-2.root");
  tree->AddFriend("gbp3", "gbp-eur-1999-2005-3.root");
  tree->AddFriend("gbp4", "gbp-eur-1999-2005-4.root");

  tree->AddFriend("jpy", "jpy-eur-1999-2005.root");
  tree->AddFriend("jpy1", "jpy-eur-1999-2005-1.root");
  tree->AddFriend("jpy2", "jpy-eur-1999-2005-2.root");
  tree->AddFriend("jpy3", "jpy-eur-1999-2005-3.root");
  tree->AddFriend("jpy4", "jpy-eur-1999-2005-4.root");

  //
  // creiamo un oggetto della classe "rete neurale"
  // (http://root.cern.ch/root/html402/TMultiLayerPerceptron.html)
  // argomenti:
  // - stringa descrittiva della rete neurale
  //   - variabili di input
  //   - numero neuroni strato nascosto
  //   - (possono esserci anche 2 strati nascosti)
  //   - variabile corrispondente all'uscita desiderata dalla rete
  // - puntatore n-upla
  // - condizione logica per selezionare campione addestramento
  // - condizione logica per selezionare campione di controllo
  //
  TMultiLayerPerceptron *mlp =
    new TMultiLayerPerceptron("usd1.q,usd2.q,usd3.q,usd4.q,"
			      "gbp1.q,gbp2.q,gbp3.q,gbp4.q,"
			      "jpy1.q,jpy2.q,jpy3.q,jpy4.q,"
			      "usd.rwi,gbp.rwi,jpy.rwi"
			      ":6:"
			      "@usd.q",
			      tree
			      // ,
			      // "Entry$%20<10 && Entry$<1000"
			      // ,
			      // "Entry$%20>=10 && Entry$<1000"
			      );

  //
  // si definiscono le liste di eventi da usare per
  // - addestramento
  // - controllo
  //
  TEventList train;
  TEventList test;
  for (int i = 0; i < tree->GetEntries(); i++) {
    if (i >= 10000) continue;
    if (i %500 < 250)
      train.Enter(i);
    else
      test.Enter(i);
  }
  train.Print();
  test.Print();
  mlp->SetTrainingDataSet(&train);
  mlp->SetTestDataSet(&test);

  // --- metodo di addesramento
  // mlp->SetLearningMethod(TMultiLayerPerceptron::kStochastic);
  // mlp->SetLearningMethod(TMultiLayerPerceptron::kBatch);
  // mlp->SetLearningMethod(TMultiLayerPerceptron::kSteepestDescent);
  // mlp->SetLearningMethod(TMultiLayerPerceptron::kRibierePolak);
  // mlp->SetLearningMethod(TMultiLayerPerceptron::kFletcherReeves);
  mlp->SetLearningMethod(TMultiLayerPerceptron::kBFGS);

  //--- quante volte iterare l'addestramento
  Int_t ntrain = 200;
  //--- inizializza i pesi delle sinapsi con valori casuali
  mlp->Randomize();

  //--- addestramento
  mlp->Train(ntrain, "text,graph,update=10");

  //
  // ora si analizzano i risultati dell'addestramento
  //

  // Use TMLPAnalyzer to see what it looks for
  TCanvas* mlpa_canvas = new TCanvas("mlpa_canvas",
				     "Network analysis");
  mlpa_canvas->SetWindowPosition(160,80);
  // mlpa_canvas->Divide(2,2);
  TMLPAnalyzer ana(mlp);
  // Initialisation
  ana.GatherInformations();
  // output to the console
  ana.CheckNetwork();
  // mlpa_canvas->cd(1);
  // shows how each variable influences the network
  // ana.DrawDInputs();
  // mlpa_canvas->cd(2);
  // shows the network structure
  mlp->Draw();
  // mlpa_canvas->cd(3);
  // draws the resulting network
  // ana.DrawNetwork(0,"type==1","type==0");
  // mlpa_canvas->cd(4);
  
  std::cout << "nnFinEst done" << endl;
}
